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In the paper we study some numerical solutions to Volterra equations which 
interpolate heat and wave equations. We present a scheme for construction of 
approximate numerical solutions for one and two spatial dimensions. Some solutions 
to the stochastic version of such equations (for one spatial dimension) are presented 
& ' as well. 
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1 Introduction 



We consider the following integrodifferential equation (Volterra type) 



fix, t) — g(x) + / a(t — s) A f(x, s) ds 



o 



t a-l 



where A is Laplacian and a(t) = — — -, T is the gamma function, 1 < a < 2, x E Mr, t > 0. 

r(a) 

The equation was considered in context of the heat conduction with memory [3,6]. 

For particular cases a = 1 and a = 2 the equation after taking the first and 
the second time derivative, becomes the heat and the wave equation, respectively. For 
1 < a < 2 the equation (JTJ interpolates the heat and the wave equations. The equation (JTJ 
was discussed extensively by Fujita [2] and Schneider & Wyss [9]. Fujita [2] has found the 
analytical form of solutions f(x,t) to (JTJ) in terms of resolvents or fundamental solutions 
S(t). 

A stochastic version of the equation ((TJ) 

f(x,t) = g{x) + f a(t-s)Af(x,s)ds + W(x,t), (2) 
Jo 

where W is some stochastic process has been studied in [4] and [5]. 

Within the resolvent approach the mild solution to (0) is given in the form: 

fix, t) = S(t)f(x, 0) + f Sit - r)dWix, t) , (3) 

Jo 

where the operator S(t) is the resolvent (fundamental solution) to the equation (JTJ), i.e. 
f(x,t) = Sit) fix, 0). The resolvent S(t) found by Fujita [2] is given by the formula 

POO 

(S(t)f)(x)= <j> ct {t,x-y)fiy)dy, t > 0, x eR, (4) 



where 

(j) a (t,x) = e" x2/4 '/v / 47rt for a = 1 and 

<f> a (t,x) = -(So(t — x) + S (t + x)) for a = 2, (5) 

(<5o(x)-Dirac's 5-function). For 1 < a < 2, the analytical form of <p a it,x) is given in 
terms of inverse Fourier transform of Mittag-Lefner function ML a iz) [2,8] and a direct 
calculation of both solutions to (JTJ) and resolvents becomes very difficult. It seems that 
obtaining approximate numerical solutions may be more practical. 

The aim of the paper is to construct: 
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• approximate numerical solutions to equations (JTJ), (j2J) (deterministic and stochastic) 
for d — 1, 

• numerical solutions to equation (P) for d = 2. 

The existing analytical solutions to (0) for d = 1 will serve as a reference to control the 
quality of the numerical approximation. 

For arbitrary 1 < a < 2, the resolvent operator S(t) for (0) does not possess a 
semigroup property. Hence, the time evolution from to t can not be divided into smaller 
steps and has to be calculated in one step. Therefore, the Galerkin method for numerical 
approach is a reasonable choice. 

The paper is organized as follows. In section 2 the Galerkin method for solving 
with one spatial dimension is presented. The numerical solutions for a = 1 and a = 2 
are compared to existing analytical ones. Examples of numerical solutions for stochastic 
cases with a simple stochastic process are presented as well. In section 3 the Galerkin 
method for two spatial dimensions is presented. Several results of numerical solutions to 
(P) for different a are shown, too. 



2 Galerkin method, case d = 1 

In Galerkin method one introduces a complete set of orthonormal functions {<pj}, j = 
1,2, ...,oo on the interval [0,t], spanning a Hilbert space H. Then the approximate 
solution is postulated as an expansion of the unknown true solution in the subspace H n 
spanned by n first basis functions {<fik}, k = 1,2, . . . ,n 

n 

f n (x,t) = ^C fc (x)0 fc (t) . (6) 



fc=l 



Inserting © into we obtain 

f* d 2 

f n (x, t) = f(x, 0) + J a(t - s)— f n (x,s)ds + e n (x,t) , (7) 
where the function e n (x,t) represents the approximation error. From © and ((Zj) we have 



e n (x,t) = f n (x,t) - f(x,0) - I a(t- s)—f n (x,s)ds 



f d 
Vc t M-/(x,0)- a(t-s) — V 

k=i Jo ax k=x 



dx 1 

c k (x)4>k(s) ds . 
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Requirement that the error function e n (x,t) has to be orthogonal to the subspace H n , 
(4>j(t),e n (x,t)) = 0, for j = 1,2, ... ,n, leads to the set of coupled differential equations 
for the coefficient functions Cj(x) 



gj(x) = cj(x) a 



k=l 



d 2 c k (x) 
dx 2 



where 



a ■, 



a(r — s)4>k(s)ds 



dr, (in general ^ a k j 



and 



gj(x) 



f( X ,0)<f>j(T)dT = f(x,0) / <f>j(T)dT. 

Jo 

Discretizing second derivative (Laplacian) one obtains as: 

1 n 

9j{Xi) = Cj(Xi) + — ^ a jk [-Cfc(Xi-i) + 2c k (xi) - c k {x i+ i) 



(9) 

(10) 
(11) 



(12) 



k=l 



with h = Xi — and j = 1, 2, . . . , n, i — 1, 2, . . . , m. 

The set (II 2j) can be written in matrix form: ^Ac = g>, where c and (7 are (N = n ■ m)- 
dimensional vectors and matrix A has a block form 



( Ci \ 

c 2 



/ [A IX ] . . . [A ln ] \ 
[An] ••• [A 2n ] 

\^ [A n l] . . . [A nn ] J 



(13) 



In (Jini) Cf = Ci(xi),Ci(x 2 ), • . . ,Ci(x m ), Gf = ^(xi), gi(x 2 ), • • . ,gi(x m ) and each block 
[Aij] is a tridiagonal matrix 



^ ^ij h 2 ®V h 2 ^ 

^ ~ 2 1 



[A, 








^ 2 1 









\ 






V 


















h 2 



1 S" 2 

2&ij Oij -\- J^dij J 



In general A is real, non-symmetric matrix (because ^ aji, see (JHIJ)). 
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2.1 Examples of numerical results for d = 1. 



Because solutions to equation $IJ are traveling wave-like functions we use free boundary 
conditions and large enough grid (precisely, f(x, t) x -+± 0O — >0 for any finite t). 



It 

or 



As initial condition we take a Gaussian distribution f(x,0) = g(x) = exp 
can represent the initial distribution of the temperature for the heat equation (a 
initial displacement of the medium for the wave equation (a = 2). 

For one spatial dimension N = n ■ m < 10 4 is usually sufficient for obtaining a 
reasonable approximate numerical solution. For such N the set of linear equations (fT2"|) 
can be solved by standard methods (e.g. LU decomposition). In fig. ^ we show numerical 
solutions to (0) for a = 1, |, |, \ and 2, at two particular time instants t = 6 and t = 12. 
The value of a in the initial condition was taken as a = 1. The reader can easily see a 
transition from a diffusion-like solution for a — 1, through intermediate cases for 1 < a < 
2, to a wave-like solution for a = 2. 

The knowledge of the analytical form of solutions for a = 1 and 2 allows us to keep 
approximation errors within a required range. To maintain the errors e(x, t) = \f&nai{x,t) — 
frmm{x,t)\ < 10~ 3 it was enough, for t — 6, to take into account a grid of m = 151 points 
in x-coordinate, covering the interval x G [—15, 15] and subspace H n with n = 8. For 
t = 12 case and the same error bounds the grid had to be increased to m = 201 points for 
the interval x G [—20, 20] and subspace H n to n = 18. Fig. 13.31 presents approximation 
errors e(x, t) for t = 6 and a = 1 and 2. 

For larger times the number of grid points and size of subspace H n has to grow in 
order to keep the same precision of numerical solutions. As the matrix A is sparse (among 
n 2 ■ m 2 elements of A at most n 2 (3m — 2) are non-zero) iterative methods for solving (fT2"|) 
become necessary. 

For stochastic equation (J2J we need some assumptions for the process W. For the first 
attempt we assumed that the process W is uniform in time, i.e. W(x, t) = CWi(x, t) (the 
constant C represents a 'strength' of the stochastic forces). Then we can approximate the 
convolution J Q * S(t — s)dW(s, x) in (J2J) in the following way: 

nt f-1 

/ Sit - s)dW(s, x) = V S(t - Si )[W(s i+lj x) - W( Si> £)], (14) 



i=0 

i 

1 ' 



where the time interval [0, t] was divided into a time grid {U = ir, i — 0, 1, ...,/}, r 
For cases a = 1 and 2, when S(i) is known analytically (see (p£J) and (JSJ) the stochastic 
convolution can be computed numerically. Fig. El compares the time evolution of solutions 
obtained numerically for a = 2 and t G [0, 6]. The top part represents the solution of the 
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deterministic equation (fTj) . the bottom one an example of a single stochastic trajectory 
(solution of the stochastic equation (J2J) with W uniform in time, C = 0.1). For more 
details and examples of numerical results, see [5]. 



3 Galerkin method, case d = 2. 

For d = 2 the equation (fT2"j) reads 

1 n 

gj(xi, yi) = Cj(xi, Vi) + a i k i~ c k( x i-ii Vi) ~ c k(%i, Vi-i) (15) 

fc=i 

+4c fc (xi, j/,) - c fc (x i+ i, y,) - c fc (xi, , (16) 

with j = 1, 2, . . . , n, % = 1, 2, . . . , m, I = 1, 2, . . . , m. 
Now, in matrix equation 

Ac = g, (17) 

c and (7 are N = n ■ m 2 -dimensional vectors, such that 
c T = ci(xi,yi), . . . ,c 1 (x 1 ,y m },c 1 (x 2 ,yi) . . . , ci(rr 2 ,y m ), • • • ,c n (x m} y m ), 
f = gifa, j/i), . . . , y m ), </i(aJ2, • • • , #1(^2, 2/m), • • • , Qn^m, Vm) 

and 





( \M ■ 


■ ■ [A ln ] 


\ 


A = 


[A21] ■ 


■■ [A 2 n] 






V [4a] ■ 


[A nn ] 


J 



(18) 



Now, every block [Aij] is the tridiagonal matrix composed of smaller blocks 



[4;] = 



/ 


(«)«■ 




(0) 


(0) 


(0) 




(0) 


\ 










(0) 


(0) 




(0) 






(0) 




(«)« 




(0) 




(0) 






(0) 


(0) 




(0) 


(«)« 


(fly 


(0) 






(0) 


(0) 


(0) 












V 


(0) 


(0) 


(0) 


(0) 






0)u 


/ 
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Blocks (ctij) are tridiagonal 



/p _|_ 4 

~h? a ij 



V 



o 
o 



ft2 a ij 
A I 4 

~h? a ij 








/i2 a «i 

























h 2 dij Sij + h 2 ciij h 2 a 



h 2 

-1 r _|_ 4 

T^T a i j Oij-rjp &i j J 



(20) 



blocks (Pij) are diagonal 



/ 



(Pij) = 



V 



o 
o 

o 
o 



h 2 dij 

































^ 









h 2 dij 



(21) 



h 2 °*i / 



and (0) are zeros, each of size m ■ m. 

Dimension of vectors c and g is N = n-m 2 . Already for n = 16, m = 200, N becomes 
large, reaching value of N = 6.4 • 10 5 unknowns and the number of matrix elements for 
the matrix A reaches N 2 = n 2 m 4 > 4 • 10 11 . 

Fortunately, the matrix A is sparse. Blocks [Aij] have at most m(3m- 2)+2(m — l)m 
non-zero elements. Then the number of non-zero elements of matrix A is at most (some 
could be 0) 

W< n 2 m(5m-4). (22) 

For n — 16, m — 200, 

W<5.12-10 7 (23) 

The size of A and its sparseness property makes using iterative methods for d > 2 
necessary. 



3.1 Conjugate and Bi-Conjugate Gradient Method 

If matrix A is symmetric and positive definite, then the problem Ax = b is equivalent to 
minimizing the function f(x) = ^x-A-x — b-x. This function is minimized when its 
gradient V/ = A ■ x — b is zero. In such case a variant of Conjugate Gradient Method is 
applicable [1, 7, 10]. 
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In our problems the matrix A is not symmetric. Therefore a more complicated Bi- 
Conjugate Gradient Method is required [1,7,10]. 

3.2 Preconditionig 

The convergent rate of iterative methods depends strongly on spectral properties of the 
matrix A. Usually matrix A is ill-conditiond. The condition number k = A max /A m i n 
is big (Aj denotes an eigenvalue). Then the convergence of iterations is usually so slow 
that accumulation of numerical errors often makes obtaining the solution impossible. The 
remedy is preconditioning. Suppose that M is a matrix that approximates A, but easier 
to invert. We can solve A ■ x = b indirectly by solving M~ l ■ A ■ x = M _1 ■ b. If 
k(M~ 1 A) <C k{A), the number of iterations is reduced significantly. 

There are several ways of choosing a preconditioner matrix M. In our case we can 
take an advantage of knowing detailed structure of matrix A (jl8Kj20|) which all elements 
are related to the elements of small matrix a ()10|) . Blocks [Aij] f)19j) and (a^) ([20)1 
are tridiagonal. We choose the preconditioner matrix M in the same block form as the 
matrix A, but leaving only diagonal blocks (<%) in |T9j) and diagonal elements 6{j + p-Oy 
in (|2()jl . All other elements of (ay) and (fly) are set equal zero. Then the matrix M 
has block form with diagonal blocks containing the same element 7^ = 5ij + mciij on 
their diagonals. Hence the matrix M _1 has the same block structure, with elements 7 _1 
on block's diagonals. The size of 7 is only n ■ n, so 7 _1 can be calculated easily by 
standard methods with the machine precision. The resulting matrix M~ X A has usually 
the condition number several orders of magnitude smaller than that of the original matrix 
A . In calculations leading to results presented below such kind of preconditionig allows 
to obtain a reasonable accuracy within 10 2 — 10 4 iterations for problems with ~ 10 5 
unknowns. 

3.3 Numerical results 

Solutions to deterministic equation for d = 1 and d = 2 differ substantially from 
each other. In fig. 0] we present the numerical solutions to for 2 spatial dimensions 
in a way convenient for comparison with fig. ^ (top) presenting solutions to (0) for 1 
spatial dimension. The initial condition for results displayed in fig. 0] is in the form 
f(x,y,0) = exp [— - ^ ] , a —2. The curves in fig. |U represent cuts of solutions f(x,y,t) 
along y — 0, i.e. f(x, 0,t) and it is clearly seen that for all given values of a the profiles 
of the solutions for d — 1 and d = 2 are different. 



8 



In fig. El two examples of the solutions f(x, y, t) at t — 6 are displayed. In the upper 
part the case a = 2 and radially symmetric initial condition is shown. In the lower part the 

case a = | with radially asymmetric initial condition (f(x,y, 0) = exp[— iiz|L] ; 

with <r 1 = 4, cr 2 = 2) is presented. 
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Figure 1: Numerical solutions to deterministic Volterra equation (P) at t — 6 (top) and 
t = 12 (bottom). Presented are cases with a = 1,|,|,|,2. Thin red line represents 
f{x,t = 0) with a = 2. 



10 




11 



X(t,x) 




Figure 3: Time evolution of the solution to the equation (1) with a = 2 for t G [0, 6], the 
deterministic solution (top) and the solution to the equation (2), a particular stochastic 
trajectory (bottom). 
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X 

Figure 4: One dimensional sections of two dimensional numerical solutions to deterministic 
Volterra equation Q at t = 6. Presented are cases with a = 1, f , |, 2 for radially symmetric 
initial conditions. 
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f(x,y,t) a=2, t=6 




f(x,y,t) a=7/4, t=6 
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Figure 5: Examples of numerical solutions to deterministic Volterra equation Q for 2 spatial 
dimensions. Presented are cases with a = 2, radially symmetric initial condition (top) at t = 6, 
and a = |, asymmetric initial condition (bottom) at t = 6. 
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